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I ' Abstract We analyze a systematic algorithm for the exact computation of the 

"t^ I current cumulants in stochastic nonequilibrium systems, recently discussed in 

the framework of full counting statistics for mesoscopic systems. This method 
is based on identifying the current cumulants from a Rayleigh-Schrodinger per- 
"jrt i turbation expansion for the generating function. Here it is derived from a simple 

C^ ' path-distribution identity and extended to the joint statistics of multiple currents. 

^ [ For a possible thermodynamical interpretation, we compare this approach to a 

'^ ■ generalized Onsager-Machlup formalism. We present calculations for a boundary 

^ ' driven Kawasaki dynamics on a one-dimensional chain, both for attractive and 

O . repulsive particle interactions. 
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1 Introduction 

r — I . We want to present and to illustrate a systematic scheme for the in principle ex- 

C^ • act computation of all possible current cumulants in Markov dynamics satisfying 

OO ! local detailed balance. The algorithm is based on an identity between current and 

^^ ■ activity fluctuations, connecting the time-antisymmetric with the time-symmetric 

Marco Baiesi 

K^ , Instituut voor Theoretische Fysica, K. U. Leuven, B-3001 Leuven, Belgium 

%_] . E-mail: marco@itf.fys.kuleuven.be 

^ ' Christian Maes 

Instituut voor Theoretische Fysica, K. U. Leuven, B-3001 Leuven, Belgium 
E-mail: christian.maes@fys.kuleuven.be 

Karel Netocny 

Institute of Physics AS CR, Prague, Czech Republic. 

E-mail: netocny@fzu.cz 



fluctuation sector as is typical for a dynamical large deviation theory in nonequi- 
librium systems. We concentrate here however on the mechanical aspect of the 
method, how it can be seen as a modified Rayleigh-Schrodinger expansion with 
specific computable expressions of the cumulants. Its relevance is therefore in re- 
liably producing also higher-order cumulants that can then be further analyzed 
for understanding the physics of some particular model. We will start that for an 
interacting particle system with boundary driven Kawasaki dynamics. 

As here we choose to emphasize the algorithm rather than its detailed numeri- 
cal implementation, we focus on relatively small systems. Yet we feel excused for 
the moment as exactly small open systems and their nonequilibrium fluctuations 
have been in the middle of attention in the last years. They are intrinsically of rel- 
evance to nanoscale-engineering and for certain cellular and molecular biological 
processes lilEI- Charge transport in nano-electromechanical systems is often de- 
scribed in terms of Markov evolutions, and is a subject of very active research 13] 
|4l|5]|6l|7]l8l. First experiments were limited to measuring the mean current or its 
variance at most, but now also third and higher order cumulants have become 
available, providing important information on quantum transitions [i9ulOI. For life 
processes, for instance in molecular motors or for ion-transport through membrane 
channels, one easily reaches energy scales as low as a few k^T II11II12II . Besides 
these cross-disciplinary aspects, the study of all these commonly called meso- 
scopic systems is important to unravel the structure of nonequilibrium statistical 
mechanics itself. 

Fluctuations cannot be ignored for small systems but rather carry signatures 
of important physics. The computational challenge in this case is not so much 
to reach large system sizes, but, for a fixed system, to obtain the fullest possible 
fluctuation patterns of the quantity of interest for long time scales. Our results 
contribute to the larger effort of organizing the computational side of the recent 
advances in nonequilibrium physics, cf. II13II141 . These theoretical results have 
often to do with fluctuation theory, as in the Jarzynski-Crooks relations 1 15,161 
or in the fluctuation theorems for the entropy production lll7lll8lfT9ll20ll21ll22ll23l . 
and going reliably beyond Gaussian characteristics in the nonequilibrium statistics 
is just a necessary but often nontrivial prerequisite. 

One of the traditional approaches to nonequilibrium solid state problems is 
the Keldysh formulation in terms of nonequilibrium Green's functions II24II25I . 
Currents of any type are obviously among the most important observables and 
their fluctuations are written in the cumulants. The basic method of the present 
paper comes from calculations within full counting statistics for small quantum 
systems ll6l l26ir7ll27ll28ll29ll30|[3Tl . We propose yet another derivation for classical 
stochastic systems based on a single path-distribution identity, which allows to 
discuss also joint current fluctuations. What follows can be seen as an adaptation to 
the framework of Markov dynamics for the computation of joint current cumulants 
in classical interacting particle systems. In Ii21ll32|[33]l34ll35ll36ll37|[38l one finds 
similar treatments. Moreover, our computational scheme aims at the same goal as 
in 09', 40, 411. but it yields exact results for small systems. The core idea of the 
method is a sufficiently simple identity, ([Hi below, that we use in combination with 
expansion techniques for eigenvalues. The novelty in our work is then as follows: 

1) We use a nonequilibrium version of the Rayleigh-Schrodinger (RS) expan- 
sion to obtain a systematic cumulant expansion for the current statistics, general- 



izing the approach of ||26| to also include joint fluctuations of different currents. 
This is particularly useful for the numerical evaluation of higher-order cumulants 
(say from third-order on) as finite-difference calculations generate more numeri- 
cal errors. One hopes that they are also within reach of experimental methods on 
real nonequilibrium systems |9 , 10|. In that case they would be invaluable tools in 
any attempt of reverse engineering. The RS expansion for (in general irreversible) 
Markov dynamics is computationally useful as every order in the expansion em- 
ploys the same basic information about the dynamics. Since the generators of 
stochastic dynamics are not symmetric matrices, their set of eigenvectors might 
be incomplete. This is taken into account in the solution of the problem, which 
involves the use of the group pseudo-inverse of stochastic matrices Ii42i . 

2) In Section |4] we add a thermodynamic interpretation of the numerical pro- 
cedure in terms of the time-symmetric sector of nonequilibrium fluctuations. This 
lines up with the recent introduction of the novel concept of traffic, which, roughly 
speaking, measures the amount of dynamical activity in the system. This activity 
counts the number of all jumps irrespectively of their direction and hence it is 
symmetric under time reversal 1.43^44^451 . It has also been considered in II37II46I . 

3) We illustrate the procedures for a boundary driven Kawasaki dynamics, for 
which an exact or analytic solution is far beyond reach, see Section |3] It is inter- 
esting to discover systematic tendencies in the role of attractive versus repulsive 
potentials for the current statistics away from equilibrium. 

The paper is organized as follows: In the next section we explain some basic 
identities that lead to the formulation of the problem as an evaluation of a certain 
eigenvalue. Section [3] gives an explicit example where the method is applied to 
a boundary driven interacting particle system. Section |4] reflects further on the 
method from the point of view of the theory of large deviations: we point out the 
role of a novel concept, that of traffic, in the interpretation of the various terms. 
The paper closes with Appendices giving details on the method and its numerical 
implementation. 



2 Method 

2.1 Current fluctuations 

We suppose a continuous time Markov jump process (X(),>o on a finite state space 
K with M elements. The dynamics is specified by all transition rates k{'r] ,^), from 
each state t] to each other ^ 7^ 77, as summarized in the generator L, which is a 
M xM matrix with elements 

Lr,r, = ~l,k{j],^) (1) 

Note that the diagonal elements equal minus the respective escape rates. We as- 
sume irreducibility in the sense that all states are reachable from any other state in 
some finite time. Hence, there is a unique stationary distribution p. We are mostly 
interested in breaking the detailed balance condition, driving the process outside 
equilibrium; see Appendix IaI for some formulation. 



The matrix L generates the stochastic evolution in the sense that 

|(/(X,)) = ((L/)(X,)> 

for all vectors / : ^ ^ M. The brackets (•) are averages both for the random (as 
yet unspecified) initial conditions as over the stochastic trajectories. The Markov 
process (X,) is a jump process in the sense that the trajectories are piecewise con- 
stant (in time) with jumps Z, = 77 — > X,+ = ^ from some state 77 to some state ^ 
at random moments t. We consider the stationary process starting from p. 

We consider each ordered pair of connected states b = (t],^) and its inverse 
is —b = (^ , 77). For a given trajectory CO = {Xt,0 <t < T), over some time-span 
T we have a microscopic current dJi,{t) = +1 when the state jumps at time t over 
the bond b, while dJi,(t) = — 1 when the state jumps over the bond —b. The time- 
integrated current 



Jo 



Jb{(0) = / dJb{t) (2) 

Jo 

thus counts the number of net jumps over the oriented bond b in the time span 

[0, T]. Note that the dependence on T in the left-hand side of (O is not explicitly 

indicated. If we look at the stationary state, we should take the expectation of (|2]l 

and divide by T to get the flux (per unit time). In the stationary state, the expected 

current over the bond b and per unit time equals 7/, = p (77)^(1],^) — p(^)fc(^,77). 

The main reason to consider all these currents like in (|2]l on the finest scale of 

transitions and the complexity of the full joint fluctuations, is to be able to move 

to arbitrary and more coarse-grained scales of description. In applications, the 

physical currents are all obtained by combinations of these currents over bonds. 

For example, an interesting current in a lattice system might count the passage of 

particles from one given site to another. In this case the current is rather of the 

form 

/b = £ 4 (3) 

beB 

where B then includes all Z? = (77 , ^ ) from a state r\ with a particle in the first site 
to a state ^ where the particle has moved to the second site (the ensemble B 
includes all bonds —b of the opposite transitions). This will in fact be our main 
example (section|3]l. 

We can formalize that: To keep the discussion as general as possible, but with 
an eye on the actual application, we consider a partition of all ordered bonds (or 
connections) ^'s consisting of sets S,B', . . . for which S, B,B', —B', . . . are mu- 
tually non-overlapping. The fully microscopic description is recovered when each 
of these sets contains exactly only one bond. 

We are interested in understanding and computing the joint fluctuations of the 
currents /g, properly rescaled as T f 0°. So for example, we want to determine the 
covariances 

C|B. = ^[(/fi^B')-W<-/B')] (4) 

in the large time T limit for B and B', which corresponds to the steady state regime. 
From now on, the bracket-expectations (• • ■ ) refer to the mathematical expectation 
in the assumed unique stationary process. Higher-order cumulants are also impor- 
tant, for example to determine the deviation from Gaussian behavior. 



In general, the computation of these cumulants like in ^ involves detailed 
information about the time-autocorrelation functions. This information is hidden 
in the spectrum of the generator L. What we will do, amounts to extracting that 
information from a systematic numerical scheme. As a further result expressions 
are obtained for these cumulants in terms of expectations of specific single-time 
observables under the invariant distribution, which allows also to see relations 
between the various cumulants and what governs them. 



2.2 General identity 

The method of computing the cumulants for the currents starts from a general 
identity ^ that relates the current fluctuations with fluctuations of occupation 
times. 

We fix a set of numbers a = (gb)- The cumulant-generating function for the 
joint fluctuations of the selected currents Jg is then given by 

gr(a) = ^log(e^^''^*> (5) 

By definition, the derivatives of gr(c7) at a = give us all possible cumulants. 
For example the second (partial) derivatives with respect to Gb, Ob' give ^. We 
therefore want to obtain an expression of gf^a) as a Taylor-expansion in the Og's, 
for T -^ +00. 

In order to do so and given the original Markov process with rates k{r\ , ^ ) we 
now construct a new Markov process with generator Lc, where the rates relative 
to bonds b= {t},^) eB and -^ = (^ , T] ) G -5 are 

£{ri,^) = k{ri,^)e'''< 

We further define the vector 

V{ri)=l^ £ fc(77,^)[e±'^«-l] (7) 

where the sign in the exponent depends on whether (t],^) = ±Z? for a selected 
bond b. 

The generating function (|5]l can be rewritten via the identity 

where the last average is over the Markov process with rates l{r\,E,), hence de- 
pending on a. 

To prove (|8]l we note that in going between the two averages (•) and (•)c7 there 
is a density e^lra)^ 

{F{(o)) = (F(o;)ee(»))^ (9) 

that is given by 



where the first sum is over all jump times t in ft) where the state changes Z, -^ X,+ , 
see for example Appendix 2 of II47I for mathematical details. As a consequence 
and via ^, 

Qi(o) = -yaBJB+ f V{X,)dt 
B Jo 

Substituting F = exp J^g CTbJb into (|9]l gives the resuh ^. 

We remark that Eq. ^ shows that the current fluctuations can be expressed in 
terms of occupation time fluctuations in a tilted path-space measure, see also in 
Section m It is not a new observation, see for instance 121113611 for very related al- 
though less explicitly stated considerations. First we continue with its exploitation 
for computational purposes. 

If one has only one set B with ag = /I 7^ 0, the current generating function 
simplifies to 

gf(A) = ilog(e^^«) (10) 

The identity ([8]l obviously remains valid, now with 

VB{rj) = [e'-l] ^ k{r^,^) + [e-'-l] £ fc(^,T]) (11) 



2.3 Feynman-Kac formula 

The right-hand side of ^ involves the single-time observable V, in contrast with 
a current being a double-time function. The V can therefore be taken as a potential 
(diagonal matrix) V in the following sense: given the matrix J^ — L^+Y, 

limllog(efc'n^,)dr)^^gmax (12) 

where e™" is the largest eigenvalue (in the sense of its real part) of ^. 

The asymptotic formula ( 112b is the limit of what is known as the Feynman-Kac 
formula. For our context, one finds a proof of it in Appendix 2 of 1,47 J . As a result, 
the current cumulants can be read from the Taylor-expansion of the eigenvalue 
e™ with explicitly known matrix 

^ = L + ^, ^ = La-L + N 

with ^ having non-zero elements only for the pairs (t] , ^ ) in some ±5 with Ob 7^ 
0. More precisely, for /? = (t] , ^ ) G B, 

^,^=fc(77,^)[.-«-l] 

Since we required that an ensemble of transitions B does not overlap with any 
other B' or —B', we can decompose the matrix ^ in a convenient sum of ma- 
trices Eb{ob) and E^b{ob), where each matrix Eb has non-zero elements equal 



to k{ri,^)e'^B only for (t],^) G S, and similarly each matrix E-g has non-zero 
elements equal to k{^ , 77 ) e^^^ only for (^ , t] ) e B. Thus, 

^ = £ [£B(aB) -£5(0) +£-B(aB) -£-b(0)] (14) 

We finally remark here that the maximum eigenvalue e™" is simple, which 
follows again from a Feynman-Kac formula saying that 

By the irreducibility assumption, the left-hand side is in fact strictly positive (for 
any T > and V), hence the right-hand side is a matrix with strictly positive 
entries. Therefore, the Perron-Frobenius theorem implies that J^ has a unique 
maximum eigenvalue. Moreover, the right and left eigenvectors of that largest 
eigenvalue of ^ have strictly positive coordinates. Exactly all the same is true for 
the generator L. 



2.4 Expansion 

From the previous discussion it is clear that ^ goes to zero with the Og's. More- 
over, there are no cross-terms containing mixed derivatives of M with respect 
to the a^'s. As we recognize the cumulants of the current distribution from the 
Taylor-coefficients in the eigenvalue c™", it is natural to write 






over the order in the Og's. Then, for each n = 1,2,... and B 

1 

ni 






This means that all odd terms n\J^g are the same matrix Eb{0) —E-b{0), and 

that all terms n!^"' with n even are equal to Eb{0) +E-b{0)- 

From the RS perturbation expansion, see Appendix IbI we obtain the following 
cumulants. It is important to note that the computation proceeds always from the 
same basic ingredients. The input consists of the stationary distribution p and the 
expression for the pseudo-inverse of L, see below. Then, all cumulants follow from 
an exact numerical calculation. More details on the algorithm are in Appendix IbI 

2.4.1 First order 

As needed, the formula for the first order cumulant corresponds to the expectation 
of the current, 

;s= Umi(7B) = (pj^i'V) (15) 

where we use the Dirac notation for left and right eigenvectors: (p | is the density 
giving the steady state occupation probabilities of the states, and 1 1 ) is a column 
vector of I's. They are the left and right eigenvectors of L, with maximal (always 
in the sense of real part) eigenvalue eo = 0. 



2.4.2 Second order 



The expression for second order gives the variance 



Cbb = /.im C^B = 2 (p I {^P -^^GJ^Pm (16) 



for the current jg over bonds with field Og, and the covariance (|4| 

Cbb' = ^lim Cj^, =- (p|(ifi''Gifi")|l> 

-(pi(^(;'Gifi")ii> 



(17) 



if B j^ B'. The matrix G is the pseudo-inverse of the matrix L in the sense of 
Drazin [42], see Appendix IaI 

2.4.3 Third and fourth cumulant 

For the higher-order cumulants we restrict to the condition of a single global cur- 
rent, as in dlOb and (lllb . In this case, we have a single ensemble B and the identity 
dHJ reduces to 

As a result, the analogue of (fT2l) is verified for the matrix J^ = L + M with ^ = 
Eb{^) — Eb{0) +E-b{X) — E-b{0). By expanding the exponential around A, —0, 
we write 

^ = £ A*^^(*) (18) 

and the cumulants are obtained from the scheme outlined in the Appendix |B] 
The third cumulant of the current Jb over bonds b E B is then 

d^) =3!(p| [^«-jfi^(i)G2^(')+^(')G^(')G^(i) 

n (19) 

-^(''G^(2)_^{2)(5_^(1) |l^, 

and the fourth cumulant is 

CW =4!(p| [^(4) -^(2)^^(2) _C^^(1)C2^(1)_(^.^)2^(1)C3^(1) 
+ if(')G^(2)c^(l)_^(l)C^(l)C^(l)C^(l) 

-if(3)Gif(i)-if(i)G^(3) 

-JB (^P) ^2^(1)^ ^(1)^2^(2) j 

+ ;Bf^(i)G^(')G2if(i)+if(i)G2^(i)G^(i) 



(20) 

Note the symmetry in the terms: when a sequence of matrices is not palindrome, 
there is also its reversed one. 



3 Example 

We consider a generalization of the symmetric exclusion process (SEP), in which, 
besides via the exclusion principle, particles are also interacting with their nearest 
neighbors at a finite reservoir temperature j3^^ . Let us consider a lattice gas on the 
sites {1 , . . . ,N}, where a configuration is an array of occupations, T] (/) = 0, 1 for 
I <i<N. The state space is thus K = {0,1}^, with M = 2^ different states. One 
can think of particles (and holes) hopping in a narrow and small (effectively ho- 
mogeneous) channel. The specific calculation below has been done for a relatively 
small system where N = S. We comment on size-dependence of the algorithm in 
Appendix O 

For the dynamics there are two modes of updating: In the bulk, a particle can 
jump to nearest neighbor sites. Then, the occupation over a nearest neighbor pair 
of sites is exchanged. For a transition 77 ^ ^ of this kind we take a rate of the form 



fc(77,^)=exp|-|(//(^)-//(77)) 



(21) 



where H is the energy function 

//(T]) = -e£77(077(' + l) (22) 

for some parameter e. Thus, only pairs of particles occupying nearest neighbor 
sites have an energetic contribution. 

At the boundaries one has the second kind of updating. At site / = 1 particles 
can be exchanged with an external reservoir having chemical potential oc/p. In 
the case of a particle leaving the system (t] (0) = 1 ^ ^ (0) = 0) the rate is given 
by 



fc(T],^)=exp|-|-|(//(^)-//(T]))' 



(23) 



while a particle enters into the system at site i = 1 with rate 

(24) 



^(77,^)=exp||-|(//(^)-//(77)) 



We focus on the time-integrated current / passing through the site i = I, which 
increases by 1 every time a particle enters there from the reservoir and decreases 
by 1 every time a particle leaves the system from there. As explained in previous 
sections, this is the sum of all microscopic currents Ji, over bonds b connecting a 
state T] with T] ( 1 ) = to another state ^ with ^ = 77 on all sites except ^ ( 1 ) = 
1. At the other boundary site i = N a similar structure may be imposed, with 
chemical potential a' /p. The time-integrated current /' however is defined there 
with the opposite convention, i.e., one adds one to J' when a particle leaves the 
system. With this convention, both currents j = J/T and / = J' /T have the same 
asymptotic value for a time span T ^ 00, as they both flow from left to right. 

The model is a boundary driven Kawasaki dynamics, reducing to boundary 
driven SEP for j3 = 0. This infinite temperature case is completely solved con- 
cerning current fluctuations in t48.49.l . If a = a', then it is easily checked that 
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the model satisfies the condition of detailed balance with respect to the grand- 
canonical distribution for energy (l22l) and chemical potential cc/p. If however 
a ^ o! then the system is driven out of equilibrium: the difference in effective 
chemical potential between reservoirs generates a particle current through the sys- 
tem. It is a fluctuating current and we study here its cumulants. For other models, 
similar questions have been addressed for example in II50II51I . Studies on the den- 
sity fluctuations for the boundary driven exclusion process are in II52II53 I. 

For simplicity, we set a' = and drive the system by varying only a and j3. 
The case a > thus corresponds to a reservoir that pushes particles from the left 
into the system, forcing a positive stationary current j. The case a < instead 
corresponds to a left reservoir that tends to remove particles. As we will see, the 
two situations are definitely not the mirror image of each other (unless /3 = 0). 

Since the product /3e is what matters in the transition rates, we simply set e = 1 
and we use the possibility /3 < for characterizing repulsive potentials. Particles 
instead attract each other for /3 > 0. Particle interactions very much complicate 
the model which is no longer analytically tractable. We use the above formalism 
to evaluate the current cumulants for different parameter values. Interestingly, in- 
teractions induce qualitatively novel behavior for the current statistics. 



3.1 Mean current 

The mean current j as a function of a and for several /3's is shown in figure[na). 
For a given j3, j increases with a, linearly around a = 0, as expected close to 
equilibrium. For each a > the mean current is maximal for a repulsive interac- 
tion (/3 < 0), see the two examples in figure|2a)). In general, the mean current j is 
not antisymmetric with respect to a, and its value in a can be very different from 
-yin-a. 

For /3 -^ — oo the problem can be mapped into the dynamics of non-interacting 
dimers "(0, 1)" and of O's. In this limit the system is somewhat like a SEP with 
I's replaced by dimers, and one thus expects a finite mean current. On the other 
hand, for j3 -^ +0°, particles stick together and it becomes more and more difficult 
for a hole (vacancy) to get in, to reach the bulk and finally to reach the other 
boundary of the channel. The hole essentially performs a random walk with an 
open left boundary before eventually reaching the system at the right boundary. 
If more than one hole enters into the system, there is a good chance that holes 
stick, further reducing the energy of the system, and their own mobility and 7 as 
well. Thus, for j3 ^ 0° we expect j -^ 0. These scenarios are qualitatively well 
confirmed in figure |2f a). 

For j3 = one has the well-studied driven SEP. In this case the current is 
antisymmetric with a, like all odd cumulants, because of a corresponding parti- 
cle/hole symmetry. 



3.2 Variance 

The second cumulant of the current distribution is its variance. For j3 = the 
variance is symmetric with a (as every other even cumulant), while for all other 
cases it displays a non-trivial dependence on a and /3, see figures [Tfb) and|2jb). 
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Fig. 1 (Color online) First four cumulants of the cuiTent distribution as a function of a, for 5 
values of the interaction "strength" /3 (see the legend). Note that for /3 = (SEP) odd cumulants 
are antisymmetric functions of the driving, while even cumulants are symmetric functions. This 
is due to a particle-hole symmetry, which is lost for interacting particles. 



For e" ^ 1 (see j3 = 4 in figure [Hb)) the variance, as the current, can reach very 
small values, confirming the scenario proposed above. 



3.3 Third and fourth cumulants 

As for the current, for j3 = (SEP) the third cumulant of the current is antisym- 
metric with a, see figurelUc). In general, however, it is a complicated function of 
a and j3, as evidenced by figure [2|c). For example, in contrast with the mean j, 
it can be a non-monotonous function of a. Similar arguments apply to the fourth 
cumulant, see figure [Hd) and figure [2|d). The third and the fourth cumulant also 
appear going to zero for e^ '^ \. 



4 Traffic 



Systematic perturbation techniques better be accompanied by a larger theoretical 
understanding. A major step in the analysis of the problem at hand that proceeds 
the numerical algorithm is contained in the simple path-distribution identity ^. 
On the left-hand side, this identity involves an average over paths (a in the origi- 
nal process, with probabilities Prob(a)). The respective probabilities in the tilted 
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Fig. 2 (Color online) First four cumulants of the current distribution as a function of /3, for 5 
values of the driving a (see the legend). Odd cumulants are identically zero when the system is 
in equilibrium (a = 0). 



space (with rates (|6]l) can be written as Probfj(ft)) = e~2(») Prob(ft)), with relative 
path-space action Q. Since on the left-hand side of identity ([8| we have a current 
generating function, a time-antisymmetric quantity is involved. On the other hand, 
on the right-hand side of (O only a potential V appears, i.e., a quantity depending 
only on the states and thus insensitive to time-reversal. Hence, the choice of the 
tilted Markov process is exactly such that the change in the time-antisymmetric 
part of the path-space action equals the appropriate sum over currents. This is 
why the exponent in the right-hand side of (l8| contains a time-symmetric function 
only. 

Such considerations are typical of the Lagrangian approach to nonequilibrium 
statistical mechanics as pioneered by Onsager and Machlup, [54|. Here however 
we are not even close to equilibrium. We thus move on a somewhat generalized 
formalism that remains quite simple for finite state space Markov processes. Nev- 
ertheless the structure of time-symmetric versus time-antisymmetric fluctuations 
is possibly important for nonequilibrium thermodynamics, if only to identify the 
relevant thermodynamic potentials, cf. 1143 1144 Il45ll . Such an identification proceeds 
via a dynamical fluctuation theory, in which we next situate the main identity ^. 

In order to rewrite (HI in another convenient form, we define occupation times 
as 



1 /-^ 
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that the path (0 = {Xf,0 < t < T) spends in state T] . Then, the exponent in the 
right-hand side of ^ equals J^ V{Xt)dt = TY^r^ V(t]);U^ (to) or, 

i^e^B'yBJB) = ^e^■Lr^v{^l)^ir,^^ (25) 

The current statistics is therefore obtained when one knows the large deviation 
rate function 7°^ for the occupation times, 

Prob^[^^«;,(77),V77]-e-""(P), T T +- 
for the modified dynamics (|6]l: 

lim ^ log {e^B ^B-'s ) = sup (;u • y - 7^ (;U ) ) (26) 

We have in mind here the application of the theory of large deviations as pioneered 
in 1551 for Markov processes. 

In that last variational expression ( 126b , the potential V also depends on a. Let 
us introduce the antisymmetric form a (77,^) = Og for (77,^) = b and a (^,77) = 
— a(T],^). Then, the change (|6]l from the original rates fc(T],^) to the new rates 
£(77 , ^ ) adds a further driving (in the spirit of local detailed balance): from (|7]l, the 
term 



(27) 



:^Y.i^^iArl:^)-^^l..k{n,i)] 



2 , 



is an expected excess traffic, defined for rates k as 

and similarly for rates (., see 1143 1144 l|45 l. The traffic expresses a time-symmetric 
kind of dynamical activity over the bond ^ = (77 , ^ ). In fact, all cumulants of the 
expansion in Section l2!4l contain the term 

^'^' ^ ' ' |tb for n even 
with expected current over bonds Z? G S 

jb = (p|[7?b(0)-7?_b(0)]|I) 
and with corresponding expected traffic 

Tb = (p|[£b(0)+£_b(0)]|I) 

For instance, the first term on the right-hand-side of ( 1161) is the stationary traffic 
Tg, while the second term can be interpreted as a zero-frequency autocorrelation 
function. 
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We stress that the traffic Tg is symmetric under the exchange 77 ^^ ^, while 
the current js is antisymmetric. In other words, the traffic adds a time-symmetric 
aspect to the evaluation of the dynamical activity. Finally, note that the following 
identity holds, 

(p 1^1 1) = XI \-'^B (coshafi - 1) + JB sinhas] 

B 
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A Markov generator and its normality 

The operator L that generates the Markov dynamics is a M x M matrix, and its spectral properties 
appear in the expansion for the cumulants (see more in the next appendix). It is important to 
realize some important changes with respect to equilibrium. For an equilibrium process with 
reversible distribution p > there is detailed balance, 

P('?)% =P{i)Hn 
Equivalently, the matrix 

^ni = Vpin]Lrj^ ■ 



is then symmetric and hence diagonalizable with a complete orthonormal set of eigenvectors. 
The matrix H is obtained from L via a similarity transformation H = Q^^LQ with here, and 
that is essential, a diagonal similarity matrix Q. In other words, we easily find a scalar product 
for which the eigenvectors of a detailed balance generator are orthonormal. All that need not be 
possible for nonequilibrium processes. 

A central notion here is that of normality: a matrix is normal if and only if it commutes 
with its adjoint if and only if it has a complete orthonormal set of eigenvectors. Detailed balance 
generators are similar with diagonal Q to normal matrices while nonequilibrium processes have 
generators that need not be similar to normal matrices at all. When such a generator is similar 
to a normal matrix, then it is diagonalizable and we can work with a bi-orthogonal family of 
left/right eigenvectors. The following example illustrates some of these points. 

Take the fully symmetric 3-state Markov process, i.e. with all rates equal to 1, and perturb 
it obtaining the generator 

-2-f + g l+f l-g 

\-f -2 + f-h I+/1 
1+g 1-/7 -2-g + h 

in the region |/|, |g|, \h\ < 1. The condition of detailed balance is satisfied on the surface f + g + 
h + fgh = 0. The nature of the spectrum depends on the sign oi D = fg + fh + gh: if D < 0, 
then the generator is diagonalizable and has real eigenvalues; if D = and at least one of the 
f.g or h is non-zero, then the matrix is not diagonalizable; if D > then it is diagonalizable 
with complex eigenvalues. In particular, all three cases occur arbitrarily close to the reference 
equilibrium f = g = h = 0. 

One consequence of the above facts directly concerns the expansion and calculation of the 
cumulants following the scheme of Appendix IB] We cannot simply rely on making use of some 
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of the standard tools of quantum mechanical calculation, like decomposition in an orthonormal 
basis. An important example concerns the calculation of the pseudo-inverse as in jl6t-lfT7t. 
When we still have a decomposition of the unity in terms of left/right eigenvectors, then the 
pseudo-inverse G can be obtained from 

G=(l-P)l(I-P)=£|4'")4)(pi°'| 

v— 1 6^ 

with (pi, ' I and \wy ) left and right eigenvectors of L with eigenvalue €{, ^ < 0, and where 

p=|i>(pi 

is the projection on the vector space of constant functions (therefore I — /* is the projection on 
the space orthogonal to them). In our general case, we employ the group inverse, a special case 
of Drazin inverse, see |42|. Its role for the computational theory of Markov processes has been 
advocated in L56J . The group inverse of L is the unique solution G of the equation 

LGL = L, GLG = G. LG = GL 

As will appear in the next section, and as visible already in {T6)-(T7} and (II9M20K that pseudo- 
inverse appears in the cumulant expansion. 

A final important difference between symmetric versus non-symmetric matrices (up to a 
diagonal similarity transformation) concerns the application of a variational principle to char- 
acterize the maximal eigenvalue. For example, in quantum mechanics one usefully employs the 
Ritz variational principle for Hamiltonians (Hermitian matrices) and for finding the ground state 
energy. We are not aware of an extension of that Ritz variational method or of a more general 
minimax principle to non-Hermitian matrices. The only variational characterization that seems 
to remain goes undirectly via the relation of the largest eigenvalue to a suitable generating func- 
tion, like in i fTlt , which itself obtains a variational expression in terms of a large deviation rate 
function, like in formula l l26t . 



B Rayleigh-Schrodinger expansion: the algorithm 

We give a review of the expansion that is used to compute the leading orders in the maximal 
eigenvalue. We refer to pages 74-81 in the book of Kato [59], for full details and for a rigorous 
treatment. 

The RS expansion finds its origins in quantum mechanical problems of time-independent 
perturbation theory |60 61 62 1. In contrast with the situation in quantum mechanics or with the 
case of detailed balance, we have in general no scalar product for which L has an orthonormal 
basis of eigenvectors. In many cases in nonequilibrium, we do have a bi-orthogonal family of M 
eigenvectors (instead of the orthonormal family in quantum mechanics) but it also happens that 
the generator is not diagonalizable and that we have no appropriate basis to express most easily 
the expansion. Fortunately, all that is not necessary and the expansion can proceed in a more 
general way. One simplifying feature is that the maximal eigenvalue that we need to compute 
is simple, as shown in Section [231 For the purpose of the present Appendix, we also make the 
simplification that only one Oij = cr ^ 0. 

The starting point is the M x M matrix L + .'% that we write in expansion 

if = L-h,'^=£cT*ifW, ^^°'^=L (28) 

i-O 

The unperturbed generator L has a resolvent r( k) = (L — k") ^ ' with Laurent series around k = 
given by 

^ ^:P+£ ic'"G"'+' (29) 

m=0 



K K 



for the projection P = 1 1) (p | on the eigenspace of eigenvalue zero, and G the pseudo-inverse in 
the sense of Drazin as we had in the previous section. 
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The resolvent for jf is 

r((J,K] = 



if-K 

defined for all k" not equal to any of the eigenvalues of Jf . It can be written as a power series in 
CT around l [29t : 

r{a,K) = r{K)+'£a"A"^{K) (30) 

n=l 

with 

rf"'(K:)= £ (-l)'V('f)-5*''''''-(K)^'''''----^'''"''-('<^) 

Vi+...+Vp=n 

where the sum is over all 1 < p < «, v, > 1. On the other hand, by Cauchy's residue theorem 

e{a) = Tr (t Kr( a, K)dK (31) 

27ti Jr 

for a circle F enclosing zero but no other eigenvalues of L. Upon substituting l lSOt into J3U we 
obtain 

e(c7) = !-Tr/ KrrrK:)^ [-i^r(K)l''d(i: (32) 

27r/ Jr 1^1 

where .^ of course depends on a. Since ^>'(k) = r{K)^, we have 

+ .^r{K)...Mr{Kf.'Mr{K) 

+ ...+.'%r(Kf ....'% r(K)Mr[K) 

Observe now that the trace and the integration commute so that l l32t becomes 

e(a) = !-Tr/ xf -^\-.'%r(K)Y dx 

^ ' Ini Jr p'lPdK-L ^ 'J 

and after integration by parts 

e(cT) = — Tr/ y -\-.'igr(K)YdK 

or 

e(c^) = -:r-Tr/log[l+,'^r()c)ldK: (33) 

27ti Jr 

Expanding the logarithm with ( I28t makes the expansion of the maximal eigenvalue 

.(0) = £aV"'=£a"— (34) 



for 

e('') = _LTr y i^lll^/^(^iV(K)...if(^"'r(K:)dK (35) 

We finally substitute the series j29t and perform the integral again with the residue theorem to 
get the result 



{n)^Y}_j±_ Y^ Trif'''')s(*'>...^(''"'5<*''' (36) 

ki+...+kp=p-[ 



p=l P V,+...+Vp=K 
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where 5'"' = —P and 5'*' = G*. The last formula can be written more explicitly obtaining the 
different orders C'"' = n\e^"^ as in Section IZ41 As an example, let us show how to compute the 
cumulant of order n = 2. The possible cases in the first sum of 06t are then p =\ and p = 2. 

For /? = 1 the second sum can only have V\ = 2 and k\ = 0, hence the contribution is 
— Tr^'^'s'"'. It is convenient to use the cyclic property of the trace operator TrAB = TrBA, 
and the definition 5'"' = —P to rewrite the term as TrP^^'. In general, given a set of left 
eigenvectors {pi\ and right eigenvectors \w(}, for the trace one has TrA = Y,t (PelM^t)- Here, 
the projection P on the 0-th eigenvectors ((p| and 1 1)) simphfies this term to (p| jf '^' 1 1). 

The only combinations of two numbers summing uptop = 2is(vi = l,V2 = l), while there 
are two choices [k] = l,k2=0) and (ki =0,k2 = I) summing up to p — 1 = 1 . The former case 
corresponds to iTrifOsflifC'sC" = iTr,if(i)G^("(-/') = jTr-P^OoifC', which 
is equal, according to previous arguments, to — 2(p|Jf'''G^'''|l). The same is true for the 

second term ^ Tr^^''^'^'^^''^''^, and their sum cancels the factor 1/2. Hence, overall one has 
the second cumulant given in Eq. l fT6t . 



C Numerical scheme 

We have shown that all that is required for the computation of cumulants, regardless of their 
order, is the information on the stationary distribution (p | and on the group inverse of the gen- 
erator, i.e. the matrix G. An efficient computation of G thus enables really making use of our 
formulas for the cumulants, like in Eg's.nst-lfTTt. jToJ, and l l20t . Given G and p, each cumulant 
is computed just by some matrix multiplications. The estimate of the group inverse of a genera- 
tor L is discussed in section 5 of |t56| and in |57|. In the computations carried out in this work, 
it turned out that 

G = P+{L-P)-^ 

was the most stable way of computing G for all parameter values. This formula derives most 
conveniently by using the properties of the so called fundamental matrix, see [581. 

However, for systems with a large number of degrees of freedom, it is rarely a good idea 
to directly invert matrices. Fortunately, it is also not necessary here. Note that any vector \z) = 
G\y) coincides with the (unique) solution of the equation L\z) = {I — P)\y) constrained by the 
condition {p\z) = 0, as immediately follows from observing that LG = GL = I — P and (p | G = 0. 
Hence, objects like GL' ' ' GL^ ' ' [ 1 ) or G~Lp''^ \ I ) can most conveniently be determined by solving 
a linear system of M equations with subsequently updated right-hand side. The number of such 
linear problems is fixed by the order of cumulants to be computed. This formulation also invites 
an application of fast iterative methods and various schemes to store sparse matrices in the 
memory, which enable to remarkably increase the system size. 

The second basic ingredient of the proposed algorithm is the computation of the stationary 
distribution (p|, for which one can choose among the available algorithms on the market. A 
possibility is to implement an Arnoldi scheme: the iteration of (p,+i | ^ (p, | {L + cl) converges 
to the eigenvector of (L + cl) with largest modulus, which coincides with (p | if the real constant 
c > is larger than the modulus of all eigenvalues of L. 

Let us finally stress that the estimates of cumulants obtained in this paper, besides having 
their own theoretical interest, have the advantage of avoiding the use of finite differences, in this 
case of eigenvalues of L^ obtained at different values of the parameters a. Like it is convenient 
to estimate the specific heat of a system from the variance of the energy distribution rather 
than from finite differences of the energy at different temperatures, we avoid the calculation of 
derivatives from finite differences, also because they usually hide dangerous dependencies on 
parameter step-sizes and the numerical instability connected with this. The latter is expected to 
be particularly problematic for cumulants of higher order. 
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